## Different choices of fixed effects (Results are estimated by lfe package version 2.8-7)
## These results are reported in Apendix Table A2-A5
rm(list=ls())
options(scipen=100,digits=10)
library(raster)
library(plyr)
library(dplyr)
library(lfe)
library(ggplot2)
library(stargazer)
load("data_final.RData")

########################
# Generate sample without post-shutdown period
data4=subset(data3, T<99)
########################
# Generate the short sample
data5=subset(data4, T>40)

##############################################################################################################################
########################
data_base1=data4
data_base2=data5

##################################################
## Year and Date FEs
all_AOD_long_C=felm(AOD~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year) + factor(T)|0|Plant_Code,
                    data = data_base1)
summary(all_AOD_long_C)

all_SO2_long_C=felm(SO2_MASS..tons.~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year) + factor(T)|0|Plant_Code,
                    data = data_base1)
summary(all_SO2_long_C)

all_NOX_long_C=felm(NOX_MASS..tons.~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year) + factor(T)|0|Plant_Code,
                    data = data_base1)
summary(all_NOX_long_C)

all_AOD_short_C=felm(AOD~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year) + factor(T)|0|Plant_Code,
                     data = data_base2)
summary(all_AOD_short_C)

all_SO2_short_C=felm(SO2_MASS..tons.~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.| factor(Year) + factor(T)|0|Plant_Code,
                     data = data_base2)
summary(all_SO2_short_C)

all_NOX_short_C=felm(NOX_MASS..tons.~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year) + factor(T)|0|Plant_Code,
                     data = data_base2)
summary(all_NOX_short_C)

stargazer(all_SO2_long_C, all_NOX_long_C,all_AOD_long_C,all_SO2_short_C, all_NOX_short_C, all_AOD_short_C, digits=3)

##################################################
## Year and Date and plant FEs


all_AOD_long_C=felm(AOD~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year) + factor(T) + factor(Plant_Code)|0|Plant_Code,
                    data = data_base1)
summary(all_AOD_long_C)

all_SO2_long_C=felm(SO2_MASS..tons.~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year) + factor(T) + factor(Plant_Code)|0|Plant_Code,
                    data = data_base1)
summary(all_SO2_long_C)

all_NOX_long_C=felm(NOX_MASS..tons.~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year) + factor(T) + factor(Plant_Code)|0|Plant_Code,
                    data = data_base1)
summary(all_NOX_long_C)

all_AOD_short_C=felm(AOD~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year) + factor(T) + factor(Plant_Code)|0|Plant_Code,
                     data = data_base2)
summary(all_AOD_short_C)

all_SO2_short_C=felm(SO2_MASS..tons.~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.| factor(Year) + factor(T) + factor(Plant_Code)|0|Plant_Code,
                     data = data_base2)
summary(all_SO2_short_C)

all_NOX_short_C=felm(NOX_MASS..tons.~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year) + factor(T) + factor(Plant_Code)|0|Plant_Code,
                     data = data_base2)
summary(all_NOX_short_C)

stargazer(all_SO2_long_C, all_NOX_long_C,all_AOD_long_C,all_SO2_short_C, all_NOX_short_C, all_AOD_short_C, digits=3)





##################################################
## Year and Date and plant*Year FEs


all_AOD_long_C=felm(AOD~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year) + factor(weekday) + factor(T) + factor(Plant_Code)|0|Plant_Code,
                    data = data_base1)
summary(all_AOD_long_C)

all_SO2_long_C=felm(SO2_MASS..tons.~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year) + factor(weekday) + factor(T) + factor(Plant_Code)|0|Plant_Code,
                    data = data_base1)
summary(all_SO2_long_C)

all_NOX_long_C=felm(NOX_MASS..tons.~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. |factor(Year) + factor(weekday) + factor(T) + factor(Plant_Code)|0|Plant_Code,
                    data = data_base1)
summary(all_NOX_long_C)

all_AOD_short_C=felm(AOD~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year) + factor(weekday) + factor(T) + factor(Plant_Code)|0|Plant_Code,
                     data = data_base2)
summary(all_AOD_short_C)

all_SO2_short_C=felm(SO2_MASS..tons.~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.| factor(Year) + factor(weekday) + factor(T) + factor(Plant_Code)|0|Plant_Code,
                     data = data_base2)
summary(all_SO2_short_C)

all_NOX_short_C=felm(NOX_MASS..tons.~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year) + factor(weekday) + factor(T) + factor(Plant_Code)|0|Plant_Code,
                     data = data_base2)
summary(all_NOX_short_C)

stargazer(all_SO2_long_C, all_NOX_long_C,all_AOD_long_C,all_SO2_short_C, all_NOX_short_C, all_AOD_short_C, digits=3)



##################################################
## Year*State and Date and weekday FEs


all_AOD_long_C=felm(AOD~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year)*factor(StateName)  + factor(weekday) + factor(T)|0|Plant_Code,
                    data = data_base1)
summary(all_AOD_long_C)

all_SO2_long_C=felm(SO2_MASS..tons.~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year)*factor(StateName) + factor(weekday) + factor(T) |0|Plant_Code,
                    data = data_base1)
summary(all_SO2_long_C)

all_NOX_long_C=felm(NOX_MASS..tons.~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. |factor(Year)*factor(StateName)  + factor(weekday) + factor(T) |0|Plant_Code,
                    data = data_base1)
summary(all_NOX_long_C)

all_AOD_short_C=felm(AOD~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year)*factor(StateName)  + factor(weekday) + factor(T) |0|Plant_Code,
                     data = data_base2)
summary(all_AOD_short_C)

all_SO2_short_C=felm(SO2_MASS..tons.~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.| factor(Year)*factor(StateName)  + factor(weekday) + factor(T) |0|Plant_Code,
                     data = data_base2)
summary(all_SO2_short_C)

all_NOX_short_C=felm(NOX_MASS..tons.~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. | factor(Year)*factor(StateName)  + factor(weekday) + factor(T) |0|Plant_Code,
                     data = data_base2)
summary(all_NOX_short_C)

stargazer(all_SO2_long_C, all_NOX_long_C,all_AOD_long_C,all_SO2_short_C, all_NOX_short_C, all_AOD_short_C, digits=3)


